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Abstract 

Beta-binomial/Poisson models have been used by many authors to model multivari- 
ate count data. Lora and Singer (Statistics in Medicine, 2008) extended such models 
to accommodate repeated multivariate count data with overdipersion in the binomial 
component. To overcome some of the limitations of that model, we consider a beta- 
bmomial/gamma-Poisson alternative that also allows for both overdispersion and differ- 
ent covariances between the Poisson counts. We obtain maximum likelihood estimates 
for the parameters using a Newton-Raphson algorithm and compare both models in a 
practical example. 

Key words: bivariate counts, longitudinal data, overdispersion, random effects, regres- 
sion models 



1 Introduction 

Beta-binomial models have been used by many authors to model binomial count data with 
different probabilities of success among units from the same group of study. Williams 
(1975) used such distributions to compare the number of fetal abnormalities of pregnant 
rat females on a chemical diet during pregnancy to a control group, both with fixed litter 
size. Gange et al. (1996) analyzed the quality of health services (classified as appropriate 
or not) during patient stay in a hospital using a similar approach. To analyze mortality 
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data in mouse litters with a fixed number of implanted fetuses, Brooks et al. (1997) used 
such models not only to allow for different probabilities of success among units from the 
same group of study, but also to consider overdispersion among them. Given that in many 
studies, the number of trials may not be fixed, Comulada and Weiss (2007) considered a 
multivariate Poisson distribution to model the number of successes and failures in a random 
number of attempts, illustrating their proposal with data from a HIV transmission study. 
Multivariate Poisson distribution have also been used to model correlated count data, as 
in Karlis and Ntzoufras (2003) who used such distribution to model the number of goals of 
two competing teams. 

In a study where the number of successes in a random number of trials was observed 
repeatedly, and therefore are possibly correlated, Lora and Singer (2008) consider multi- 
variate beta-binomial/Poisson models. In their proposal, the beta-binomial component also 
accounts for overdispersion across units with the same levels of covariates. The multivari- 
ate Poisson component accommodates both the random number of trials and the repeated 
measures nature of the data. The effect of possible covariates is taken into account via 
the regression approach suggested by Ho and Singer (1997, 2001). Their model, however, 
requires a constant covariance term between the repeated number of trials and does not 
allow for overdispersion in these counts. Since, as suggested by Cox (1983), the precision 
of parameter estimates may be seriously affected when overdispersion is not accounted for 
in the models considered for analysis, we propose a beta-binomial/gamma-Poisson model 
that not only incorporates such characteristics but is also easier to implement computa- 
tionally. The model, along with maximum likelihood methods for estimation and testing 
purposes are presented in Section 2. An illustration using data previously analyzed by Lora 
and Singer (2008) is presented in Section 3. A brief discussion and suggestions for future 
research are outlined in Section 4. 

2 The beta-binomial/gamma-Poisson model for repeated mea- 
surements 

We denote the vector of responses for the g-th sample unit (g = 1, . . . , M) by 
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with X g h corresponding to the number of successes in N g h trials performed under the h-th 
(h = 1, . . . ,p) observation condition. We assume that for all g and h, 

Xgh | Nghj-Kgh follow independent binomial ( N g h, -K g h) distributions (1) 
TTg h follow independent Beta(n(z^ gh )/9(zg gh ), [1 - fi(z^ gh )]/9(z egh )) distributions^) 
Ngh | T g follow independent Poisson(A(zA s /i)r s ) distributions (3) 
T g follow independent gamma(a(z a9 )/i5(zj 9 ), l/5(zs g )) distributions (4) 

where z^gh, zg g h, z\ g h, z ag and z$ g are vectors of fixed covariates. 

According to (1) and (2), the success probabilities may be different across units, but 
they are generated by beta distributions that may depend on covariates. In (3) and (4), 
we follow Nelson (1985) to specify that the numbers of trials may also be different across 
units, but are generated by gamma distributions that may also depend on covariates. 

The parametrizations (0 < fj, < 1, 9 > 0) adopted in (2) and (a > 0, 5 > 0) adopted 
in (4) are used to facilitate maximum likelihood estimation, as suggested by Gange et 
al. (1996); their relation to the usual beta(a, b) parametrization, as in Johnson and Kotz 
(1970), and the usual gamma(c, d) parametrization, as in Mood et al. (1974), is given by 

/j, = -, 9 = -, a = — and 8 = —. 

a + b a + b a a 

The first and second order central moments of T g in (4) are 

E(t 9 ) = a(z ag ) (5) 
Var(T g ) = a{z ag )5{z Sg ) (6) 

From (3) and (4), the first and second order central moments of the number of trials are 

E(N gh ) = \{z Xgh )a{z ag ) (7) 
Var(N g h) = X(z Xg h)a(z ag ){l + \(z\ gh )5(z Sg )} (8) 
Cov(N gh , N gh ,) =X(z Xgh )X(z )*M (9) 

for all g,h,h', h ^ h! . Similarly, the first and second order central moments of ix g h in (2) 
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E(TT gh ) = M(z M gfc) (10) 
Var(ir gh ) = p,{z^ gh )[l - ii(z^ gh )]9(z d gh)[l + 9{z egh )}~ 1 (11) 
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Also, from (1) and (2), we may conclude that, for all g and h, 

Xgh | N gh ~ beta - binomial [N gh , /J.(z^ gh ), 9(zg gh )] 

with 

E{X gh ) = fi(z ligh )X(zxgh)a{z ag ) (12) 

Var(Xg h ) = ^(z wfc )[l - /i(z wft )] - A 2 (z Agfe )a(z a g)[a(z ag ) + <J(z 5p )] 

+M z /^)M z Ag/Ja( z a ff )[l + v{ z Mh)K z \gh)d(zSg)} (13) 
Cov(X gh ,X gh >) = /j,(z llgh )ij l (z fJigh/ )X(z Xgh )X(z Xgh/ )a(z ag )8(zsg) (14) 

for all g,h,h' , h ^ h' . The covariance between the numbers of successes and trials is 

Cov(X gh ,N gh ) = fl(z^g h )\(z X gh)tt(z a g){l + \(z X gh)5(z6g)}. (15) 

The parameters 6(zg g h) govern both the variability of the success probabilities and the 
overdispersion of the number of successes, that may also depend on the parameter 6(zg g ). 
When 9{zQgh) and 5(zs g ) are equal to zero, there is no overdispersion for the number of 
successes. The parameters S(zg g ) are also related to the variability and overdispersion of 
the number of trials and to the covariance between the numbers of trials and numbers of 
successes. When 5(zs g ) = 0, the repeated counts are independent. 

To investigate the effects of covariates, we adopt log-linear models of the form 

" (Z " a ' J = i + »pK, 9 A) (16) 

0(z egh ) = exp(z' ggh (3 g ) (17) 

X(z Xgh ) = exp(z' Xgh f3 x ) (18) 

a(z ag ) = exp{z' ag (3 a ) (19) 

5{z Sg ) = exp(z' Sg P 5 ) (20) 

where /3„, (3 9 , (3 X , (3 a and (3g are vectors of parameters to be estimated. 

From (1), (2), (3) and (4) it follows that the joint probability mass function for the 
number of trials and successes for the 5-th unit is 

v 

P(X gl ,N gl ,...,X gp ,N gp ) = Yl P(X gh I N gh )P(N gl , N gp ) 

h=l 
P 



r 1 poo f 

[J P(X gh I N gh ) / J] P(N gh I r g )f(r g )dr g 
h=i \ Jo h=i 
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with / denoting the density of (4). Since the logarithm of the likelihood is given by 



M p M 

= l °9 p (X 9 h I N gh , 0^, fa) + logP(N gl ,..., N gp \p x ,P a ,p 5 ), 

9=1 h=l 9=1 

the parameters of the beta-binomial distribution (/3 ,/3 e ) can be estimated separately from 
those of the gamma- Poisson distribution (/3 A , (3 a , P s ). 

The beta-binomial probability mass function can be written as 



P(X gh = Xg h | Ng h = UghlPpPe) 



n gh 

\ x gh, 



1 



O(Z0g h ) 



1 



O(Z0g h ) 



+ n 



-l 



x r 



x r 



{e(z egh ) +Xgh 



\0{*0gh)) 



1 ~ /^(z M 9fe) 
0{z8gh) 



-\- Tlgh x gh 



r 



1 - K z i*gh] 

0(z9gh) 



-1 



ngH ) II [l + vfifaghT 1 II H^gh)+v9(z egh )} 
\ x ghJ u =o 



v=0 



[1 - K z /igh) + w6(z 



(21) 



w=0 



where T(r) = / °° t r ~ 1 e~ t dt. The expressions involving ratios between two gamma functions 
(presented within brackets) make sense when n g h ^ (in the first ratio), x g h ^ (in the 
second ratio) and x g h ^ n g h (in the third ratio). When these conditions are not satisfied, 
the ratios between the gamma functions may be set equal to one, and do not affect the 
conditional probability of X g h given Ngh, (3^, Pg. 

The kernel of the beta-binomial log-likelihood function is 



M p 

EE 

9=1/1=1 



L gh- 



log[fi(z^g h ) + v9(zg gh )] + 

v=0 

X] log[l - m( z m9/i) + w9(z egh )] - ^ log[l + u6(z dgh )] 



w=0 



u=0 



(22) 



and we may use maximum likelihood methods adopting a Newton-Raphson iterative process 
to estimate P^ and Pg. The first and second derivatives of (22) are shown in Lora and Singer 
(2008). Method of moments estimates based on the beta-binomial distribution may be used 
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as initial values for ^{z^gh) and 9(zg g h), as suggested by Griffiths (1973). Likelihood ratio 
tests may be employed for model reduction purposes, i.e., for constructing a parsimonious 
model that captures the explainable variability in the data. For example, to verify if the 
q-parameter vector (3* is null, the test statistics LR = 2(L — L*), with L* indicating the 
log-likelihood under Hq and L, this logarithm under the alternative hypothesis may be 
employed. Asymptotically, LR follows a chi-squared distribution with q degrees of freedom 
under the null hypothesis. 

The probability function for the repeated number of trials based in (3) and (4) is 



P(iV, 



si 



n 



si) 



.Ngp = n gp \(3 x ,(3 a ,f3 s ) 



n 



h=l 
P 



.h=l 
P 

n 

h=l 



n gh \ 

E A ( Z ^) + ^ 



1 



a{z ag )/&(z Sg ) 



T, p h=1 n gh +a(z ag )/8(z Sg ) 



J2righ + 



Jl=l 




n gh 

V 



^ P h = i n 9h- 



[a(z ag ) + u5(z Sg )} 



h=l 



Xgh) 



u=0 



+ 1 



(23) 



In (23), the simplifications for the rations between two gamma functions make sense when 
a n g h / 0. When this condition is not satisfied, the ratio is also set equal to one, and 
it does not affect the probability value. 

The kernel of the gamma-Poisson log-likelihood function is 

M ( p Z p h=1 n gh -1 

L(P\,P a ,Ps) = E { ^2l n 9hlogX(z Xgh )] + log[a(z ag ) + u5(z 5g )] 



9=1 {h=l 

a.{z ag ) 



u=0 



h=\ 



log 



6 ( z Sg) E A M + 1 



\h=l 



(24) 



and we adopt the same methods used with the beta-binomial model to estimate (3^ , j3 a and 
j3g. The first and second derivatives of (24) are shown at the Appendix. Method of moments 
estimates may be used used as the initial values for \ g h{z\), a g (z a ) and 8(z§) here, too. 
Likelihood ratio tests may be employed for model reduction purpose, along similar lines as 
those considered for the beta-binomial model. 

Both iterative processes are implemented in the R software and the corresponding code 



can be downloaded from http://www.ime.usp.br/~jmsinger 
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3 Data analysis 



To compare the beta-binomial/gamma- Poisson to the multivariate beta-binomial/Poisson 
model, we consider the same data presented in Lora and Singer (2008) from a study con- 
ducted at the Learning Laboratory of the Department of Physiotherapy, Phonotherapy and 
Occupational Therapy of the University of Sao Paulo, Brazil, to evaluate the performance 
of some motor activities of Parkinson's disease patients. For the sake of completeness, we 
repeat the description of the study here. Twenty five patients with confirmed clinical di- 
agnosis of Parkinson's disease and twenty one normal (without any preceding neurologic 
alterations) subjects repeated two sequences of specified opposed finger movements (touch- 
ing one of the other four fingers with the thumb) during one minute periods, with both 
hands. This was done both before and after a four-week experimental period in which 
only one of the sequences was trained (active sequence) with one of the hands; the other 
sequence was not trained (control sequence). Half of the subjects in each group trained the 
preferred hand (right for the right-handed and left for the left-handed in the normal group 
or the less affected by the disease in the experimental group) and the other half trained the 
non-preferred hand. Information on the number of attempted and successful trials were 
recorded with a special device attached to a computer. 

Six subgroups may be characterized by the combination of disease stage (normal, ini- 
tial or advanced) and use of the preferred hand (yes or no). The repeated measures are 
characterized by the cross-classification of the levels of sequence (control or active) and 
evaluation session (baseline or final). The specific objective of the study was to evaluate 
whether training is associated with increases in the expected number of attempted trials 
per minute (agility) and/or on the probability of successful trials (ability). Note that the 
treatment could improve agility without improving ability, so an evaluation of its effect on 
both characteristics is important. 

The means and variances of the number of attempted and successful trials at the baseline 
and final evaluations with the active and control sequences for patients at the different dis- 
ease stages using the preferred or non-preferred hands are presented in Table 1. Variances, 
instead of standard deviations, are displayed to facilitate identification of over dispersion in 
the sense referred by Nelder and McCullagh (1989), i.e., cases where variances are greater 
than expected under Poisson or binomial distributions. Overdispersion in the number of 
attempts, under a Poisson distribution is clearly identified by comparing the observed mean 
and variance; for the number of successes, on the other hand, it is necessary to compare the 
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observed and expected variances under the binomial distribution (np(l —p))- For example, 
considering normal subjects performing the active sequence at the baseline session using the 
preferred hand, the expected variance under the binomial model is 1.4, while the observed 
variance is 49.0, highlighting the overdispersion for these counts too. 

Correlation coefficients for the within-subject responses for the normal patients using 
the preferred hand are displayed in Table 2. For this subgroup, only 3 out of the 28 observed 
correlations are smaller than 0.60; this suggests that the counts are probably related and it 
is sensible to use a model that can accommodate this relationship. The correlation patterns 
for the other subgroups are similar and are not presented. 

The analysis strategy consisted in fitting initial models of the form (16)-(20) with all 
main effects and first order interactions, and trying to reduce them by sequentially elimi- 
nating the non-significant terms. The parameters are indexed by disease stage (0=normal, 
l=initial, 2=advanced), intervention hand (P=preferred, N=non-preferred), evaluation ses- 
sion (B=baseline, F=final) and sequence (C=control, A=active). We adopted a reference 
cell parameterization with the reference cell corresponding to the normal group (0), per- 
forming the active sequence (A) with the preferred hand (P) at the baseline evaluation 
(B). 

3.1 Modelling the expected probability and dispersion of successful at- 
tempts 

For both beta-binomial/gamma- Poisson and multivariate beta-binomial/Poisson models, 
the parameters of the beta-binomial components can be estimated separately from those 
of the gamma-Poisson or the multivariate Poisson distributions. Therefore, modelling the 
expected probabilities and dispersion parameters of the successful attempts is exactly the 
same as in Lora and Singer (2008) and it is not shown here; we present only the estimates 
and standard errors computed under the final beta-binomial model (Table 3) for comparison 
with the results obtained under the beta-binomial/gamma- Poisson model. Under this final 
model, estimates of the expected probabilities of successful attempts [Efagh) = fi(z w ),)] and 
dispersion parameters 6(zg g f l ) (that govern the variability of the probabilities of successful 
attempts), along with their standard errors, are presented in Table 4. 

The results suggest no evidence of difference between the expected probabilities of 
successful attempts for patients using preferred or non-preferred hand {P^n = 0), neither 
for active nor for control sequences in the baseline session (fiaC = 0). Patients in the normal 
group or with the disease in initial stage have similar expected probabilities of successful 
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Table 1: Mean and variance (within parentheses) of the number of attempted and successful 
trials. 



Disease 


Evaluation 


Intervention 


Sequence 


Successes 


Attempts 


stage 


session 


hand 








Normal 


Baseline 


Preferred 


Control 


17.1 (49.0) 


18.6 (46,2) 


Normal 


Baseline 


Preferred 


Active 


17.1 (72.3) 


17.9 (79.2) 


Normal 


Baseline 


Non-preferred 


Control 


18.1 (27.0) 


20.9 (47.6) 


Normal 


Baseline 


Non-preferred 


Active 


17.1 (37.2) 


19.5 (53.3) 


Normal 


Final 


Preferred 


Control 


20.9 (90.3) 


26.1 (44.9) 


Normal 


Final 


Preferred 


Active 


32.7 (139.2) 


33.1 (132.3) 


Normal 


Final 


Non-preferred 


Control 


24.2 (25.0) 


28.6 (38.4) 


Normal 


Final 


Non-preferred 


Active 


32.8 (74.0) 


34.4 (72.3) 


Initial 


Baseline 


Preferred 


Control 


13.7 (24.0) 


16.3 (44.9) 


Initial 


Baseline 


Preferred 


Active 


12.0 (23.0) 


13.5 (23.0) 


Initial 


Baseline 


Non-preferred 


Control 


12.0 (17.6) 


14.6 (9.0) 


Initial 


Baseline 


Non-preferred 


Active 


10.7 (20.3) 


13.6 (10.9) 


Initial 


Final 


Preferred 


Control 


13.2 (30.3) 


16.8 (43.6) 


Initial 


Final 


Preferred 


Active 


20.2 (9.6) 


21.8 (2.9) 


Initial 


Final 


Non-preferred 


Control 


15.3 (112.4) 


20.3 (116.6) 


Initial 


Final 


Non-preferred 


Active 


20.1 (33.6) 


20.4 (39.7) 


Advanced 


Baseline 


Preferred 


Control 


4.8 (22.1) 


7.1 (11.6) 


Advanced 


Baseline 


Preferred 


Active 


4.6 (11.6) 


7.9 (14.4) 


Advanced 


Baseline 


Non-preferred 


Control 


8.3 (72.3) 


12.5 (15.2) 


Advanced 


Baseline 


Non-preferred 


Active 


13.5 (92.2) 


15.5 (57.8) 


Advanced 


Final 


Preferred 


Control 


7.4 (75.7) 


11.9 (67.2) 


Advanced 


Final 


Preferred 


Active 


13.5 (90.3) 


14.9 (77.4) 


Advanced 


Final 


Non-preferred 


Control 


5.8 (31.4) 


12.8 (12.3) 


Advanced 


Final 


Non-preferred 


Active 


22.5 (75.7) 


23.8 (75.7) 
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Table 2: Correlation coefficients for the within-subject responses for the normal subjects 
using the preferred hand 











Baseline 


session 






Final 


session 








Active seq. 


Control seq. 


Active seq. 


Control seq. 








Sue. 


Att. 


Sue. 


Att. 


Sue. 


Att. 


Qnn Aff 

oUC. Att. 


Baseline 


Active 


Sue. 


1 














session 


seq. 


Att. 


0.99 


1 














Control 


Sue. 


0.85 


0.84 


1 












seq. 


Att. 


0.78 


0.80 


0.96 


1 








Final 


Active 


Sue. 


0.76 


0.76 


0.61 


0.61 


1 






session 


seq. 


Att. 


0.74 


0.74 


0.61 


0.63 


0.99 


1 






Control 


Sue. 


0.53 


0.49 


0.59 


0.63 


0.60 


0.61 


1 




seq. 


Att. 


0.81 


0.82 


0.70 


0.69 


0.93 


0.92 


0.50 1 



Codes: Sue. —Successes, Att.=Attempts and seq.=sequence 



attempts (/3^i = 0) , but those with the disease in an advanced stage have smaller expected 
probabilities of successful attempts ifi^a, < 0). Moreover, an intervention effect is detected 
since the expected probabilities of successful attempts in the final session are greater than 
those for the baseline session (/3uF > 0). These values are smaller for the control sequence 
than for the active sequence (P^f + P/i(f*c) < 0) suggesting that training is effective with 
respect to ability. 

We may also infer that there is no difference between the expected dispersion parameter 
for subjects performing the active and control sequences (Pec = 0). F° r the normal subjects, 
the expected dispersion parameters are the same (Pec, PeN, PeF=ty, except in the final 
evaluation using the non-preferred hand, for which the expected value is smaller than the 
others (Pe(F*N) < 0). For patients in initial stage of the disease, the expected dispersion 
parameters are smaller than for those in the normal group (Pe\ < 0); however, they change 
for each combination of session and intervention hand (Pe(i*F)j Pe(i*N)-> Pe(F*N) 7^ 0). 
Finally, for patients in the advanced stage of the disease, the expected dispersion parameter 
is larger than for those in the normal group (Pei > 0), but this changes for the final session 
when the non-preferred hand is used (Pe(F*N) 7^ 0). 
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Table 3: Parameter estimates and standard errors under the final beta-binomial model 









Standard 


Parameter 


Related to 


Estimate 


error 




Normal group, preferred hand, 
baseline session and active sequence 


1.86 


0.15 




Effect of advanced stage 


-1.35 


0.25 




Effect of final session 


1.38 


0.30 




Effect of final session and control sequence 


-1.79 


0.30 


f3 S 


Normal group, preferred hand, 
baseline session and active sequence 


-1.07 


0.27 


Pei 


Effect of initial stage 


-2.98 


1.05 


1392 


Effect of advanced stage 


1.31 


0.37 


Pe{i*F) 


Effect of disease in initial stage and final session 


1.66 


0.82 




Effect of initial stage and non-preferred hand 


2.78 


0.91 


f3e(F*N) 


Effect of final session and non-preferred hand 


-1.49 


0.44 



3.2 Modelling the expected number of attempts 

The initial model parameter vector, with all main effects and first order interactions is 
P = (Px^P^Ps) wnere 

Px = Wxo, Pxi, Px2, PxN, PxF, PxC, 

Px{l*F) ; Px(l*N) i Px(l*C) ) Px(2*F) ) Px(2*N) > fix(2*C) > Px(F*N) > fix{F*C) ) Px(N*C) ) 
Pm = (AnOj /^ml, Pm2i PmN, Pm(l*N) ; (3 m (2*N)) 

with m = a, 5. We may interpret (3\o as the logarithm of A for normal individuals, using 
the preferred hand, performing the active sequence at the final evaluation; (3\n corresponds 
to the variation in the logarithm of A due to the effect of the non-preferred hand compared 
to the preferred one; /^(l*^) corresponds to an additional variation in the logarithm of A 
due to the interaction between the initial stage of the disease (1) and the use of the non- 
preferred hand (N). The elements of the vector /3 A related to different evaluation sessions 
(represented by F and C) allow for different number of attempts in these different evaluation 
sessions. On the other hand, a(z ag ) and 5(zg g ) do not vary in different evaluation sessions; 
therefore the vectors j3 a and (3§ do not have elements to distinguish between sessions, but 
have elements to compare subgroups. 

As noticed in Lora and Singer (2008) for the beta-binomial model, the iterative process 
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Table 4: Estimates of expected probabilities of successful attempts, dispersion parameters 
and standard errors under the final beta-binomial model 



Disease 


Evaluation Intervention 




Expected 


Standard 


stage 


session 


hand 


Sequence 


value 


error 


Expected probabilities of successful attempts 


Normal or initial 


Baseline 


Either 


Either 


0.87 


0.02 


Normal or initial 


Final 


Either 


Control 


0.81 


0.03 


Normal or initial 


Final 


Either 


Active 


0.96 


0.01 


Advanced 


Baseline 


Either 


Either 


0.62 


0.06 


Advanced 


Final 


Either 


Control 


0.52 


0.06 


Advanced 


Final 


Either 


Active 


0.87 


0.04 


Dispersion parameters 


Normal 


Baseline 


Either 


Either 


0.34 


0.09 


Normal 


Final 


Preferred 


Either 


0.34 


0.09 


Normal 


Final 


Non-preferred 


Either 


0.08 


0.03 


Initial 


Baseline 


Preferred 


Either 


0.02 


0.02 


Initial 


Baseline 


Non-preferred 


Either 


0.28 


0.14 


Initial 


Final 


Preferred 


Either 


0.09 


0.06 


Initial 


Final 


Non-preferred 


Either 


0.33 


0.19 


Advanced 


Baseline 


Either 


Either 


1.27 


0.37 


Advanced 


Final 


Preferred 


Either 


1.27 


0.37 


Advanced 


Final 


Non-preferred 


Either 


0.29 


0.13 



was very sensitive to initial values, specially for the interactions. To overcome this difficulty, 
we started with a simpler model containing only the main effects and used the resulting 
estimates as initial values for fitting other models, obtained by including the interactions one 
by one. The estimates of the interaction parameters obtained in this preliminary process 
were used as the initial values in our modelling strategy. 

The non-significant interactions were identified and their simultaneous elimination from 
the initial model was supported (p = 0.211) via a test of the hypothesis 

Hq : P\(1*F) j P\(l*N) > P\(1*C) ! /?A(2*F) , P\(2*N) , P\(2*C) j P\(F*N) j P\(N*C) > 
Pa(l*N) > Pa(2*N) j Ps{l*N) > Ps(2*N) = 

Under the resulting reduced model, the non-significant main effects were identified; their 
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simultaneous elimination was corroborated {p = 0.493) via a test of the hypothesis 

Hq : P\N, Pxd Pal, Pa2, PaN, Psi, Ps2, PsN = 0. 

We considered other hypotheses where some of these parameters are equal to zero and they 
were all rejected {p < 0.150). Goodness of fit of the resulting reduced model was confirmed 
by a likelihood ratio test in which it was compared to the initial model {p = 0.289). 

For this final model, the corresponding parameter estimates along with their standard 
errors are presented in Table 5. Based on this, we estimated expected values for X(z\ g h); 
the results are presented in Table 6. Additionally, since only the parameters p a o and Pso 
were included at the final model, we have a(z ag ) = 3.67, with standard error of 0.18, and 
5(zs g ) = 0.27, with standard error of 0.07, for all disease stages and both hands. The 
non-zero estimate of 5 suggests that the total attempts are overdispersed and that the 
correlations among the counts across the different instants of evaluation are non-null. 

Table 5: Parameter estimates and standard errors for the final gamma-Poisson model 



Parameter 


Related to 


Estimate 


Standard error 


Pxo 


Normal group, preferred hand, 
initial evaluation and active sequence 


1.68 


0.03 


Pxi 


Effect of initial stage 


-0.38 


0.05 


Px2 


Effect of advanced stage 


-0.71 


0.05 


PxF 


Effect of final evaluation 


0.52 


0.04 


Px(F*C) 


Effect of final evaluation and control sequence 


-0.22 


0.05 


PaO 


Normal group, preferred hand 


1.30 


0.05 


Pso 


Normal group, preferred hand 


-1.32 


0.25 



We may conclude that individuals in the initial stage of the disease have smaller expected 
number of attempts than normal ones, and for individuals in the advanced stage this value 
is even smaller (P\ 2 < P\i < and p a \ = p a 2 = 0). There is no evidence of difference 
between the expected number of attempts for participants using preferred or non-preferred 
hands (P\n = and p a ^ = 0), neither for active nor for control sequences in the baseline 
session (P\c = 0). The results suggest that the training is also effective with respect to 
agility, since the expected number of attempts under the final evaluation is bigger than 
at the initial one (P\f > 0). Moreover, for the control sequence, the expected number of 
attempts is larger at the final evaluation compared with the initial one (P\f + P\(f*c) > 0)i 
however, considering only the final evaluation, the expected number of attempts is larger 
for the active sequences than for the control ones (P\(f*c) < 0)- 
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Table 6: Estimates of expected values of X(z\ g h) 



Disease 


Evaluation 


Intervention 


Sequence 


Expected 


Standard 


stage 


session 


hand 




value 


error 


Normal 


Baseline 


Either 


Either 


5.4 


0.2 


Normal 


Final 


Either 


Control 


7.2 


0.3 


Normal 


Final 


Either 


Active 


9.0 


0.4 


Initial 


Baseline 


Either 


Either 


3.7 


0.2 


Initial 


Final 


Either 


Control 


5.0 


0.4 


Initial 


Final 


Either 


Active 


6.2 


0.3 


Advanced 


Baseline 


Either 


Either 


2.3 


0.1 


Advanced 


Final 


Either 


Control 


3.6 


0.2 


Advanced 


Final 


Either 


Active 


4.4 


0.3 



Table 7 contains estimates of the expected successful and total attempts along with 
the respective standard errors. In Table 8 we present estimates (with respective standard 
errors) of the elements of the covariance matrix for normal subjects using the preferred 
hand. Covariance patterns for the other subgroups are similar and are not included. 

Table 7: Estimates and standard errors (within parentheses) for the expected number of 
successful and total attempts under the final beta-binomial/gamma- Poisson model 



Disease 


Evaluation 


Intervention 


Sequence 


Successful 


Total 


stage 


session 


hand 




attempts 


attempts 


Normal 


Baseline 


Either 


Either 


17.2 (1.0) 


19.8 (1.1) 


Normal 


Final 


Either 


Control 


21.4 (0.8) 


26.4 (0.1) 


Normal 


Final 


Either 


Active 


31.7 (1.8) 


33.0 (1.8) 


Initial 


Baseline 


Either 


Either 


11.8 (0.7) 


13.6 (0.8) 


Initial 


Final 


Either 


Control 


14.9 (1.1) 


18.4 (1.2) 


Initial 


Final 


Either 


Active 


21.9 (1.4) 


22.8 (1.4) 


Advanced 


Baseline 


Either 


Either 


5.2 (0.6) 


8.4 (0.6) 


Advanced 


Final 


Either 


Control 


6.9 (0.9) 


13.2 (0.9) 


Advanced 


Final 


Either 


Active 


14.0 (1.2) 


16.1 (1.1) 



4 Discussion 

The proposed beta-binomial/gamma-Poisson model is more general than the multivariate 
beta-binomial/Poisson model considered in Lora and Singer (2008) because it allows for 
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Table 8: Estimates and standard errors (within parentheses) for the expected covariance 

matrix for normal subjects using the preferred hand 

Baseline session Final session 

Active seq. Control seq. Active seq. Control seq. 









Sue. 


Att. 


Sue. 


Att. 


Sue. 


Att. 


Sue. 


Att. 


Baseline 


Active 


Sue. 


51.2 
















session 


seq. 




(7.2) 




















Att. 


42.2 


48.5 




















(11.4) 


(13.0) 
















Control 


Sue. 


21.5 





51.2 














seq. 




(5.8) 




(7.2) 
















Att. 





28.7 


42.4 


48.5 


















(7.8) 


(11.4) 


(13.0) 










Final 


Active 


Sue. 


40.3 





40.3 





118.9 








session 


seq. 




(10.8) 




(10.8) 




(22.2) 












Att. 





48.4 





48.4 


110.5 


115.1 














(13.0) 




(13.0) 


(30.0) 


(31.2) 








Control 


Sue. 


27.3 





27.3 





52.3 





87.4 






seq. 




(7.3) 




(7.3) 




(13.8) 




(12.9) 








Att. 





39.0 





39.0 





65.8 


64.9 


80.1 










(10.5) 




(10.5) 




(17.7) 


(17.8) 


(21.8) 



Codes: Suc.=Successes, Att.=Attempts and seq.=sequence 



different covariances between the number of attempts in different evaluation sessions and 
considers a possible over dispersion of the total attempts. Moreover, the gamma-Poisson 
component of the model is computationally much easier to use for comparisons among the 
numbers of attempts in different evaluation sessions. 

While in the multivariate beta-binomial/Poisson model, the multivariate Poisson com- 
ponent requires a different set of parameters for each evaluation session, in the beta- 
binomial/gamma- Poisson model, the gamma-Poisson component includes a single set of 
parameters for all evaluation sessions. To compare the expected number of attempts under 
different conditions using the former, it is necessary rewrite the model and to derive ad 
hoc estimating equations while under the latter, it suffices to eliminate the corresponding 
regression parameter and to obtain new parameter estimates using the same estimating 
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equations. For the analyzed data, for example, the comparison between the control and ac- 
tive sequence during the baseline evaluation using the beta-binomial/gamma-Poisson model 
is done by testing if the parameter /3\C is null. On the other hand, under the multivari- 
ate beta-binomial/Poisson approach, the total number of trials is modelled with a specific 
vector of parameters for each instant of observation; for the data in the example, they are: 
baseline evaluation performing active sequence, baseline evaluation performing the control 
sequence, final evaluation performing the active sequence and final evaluation performing 
the control sequence. To compare the control and active sequences during the baseline 
session we should rewrite the model using only three parameters: baseline evaluation (the 
same for active and control sequences), final evaluation performing active sequence and 
final evaluation performing control sequence. 

The average of the absolute differences between the sample means of the number of 
successful and total attempts and the respective expected values under this final model 
(Table 7) is 1.7. The same average based on the multivariate beta-binomial/Poisson model 
is 0.9. Furthermore, the average of the absolute differences between the observed and 
estimated covariances using the multivariate beta-binomial/Poisson model is 21.5 while 
it is 19.1 if we use the beta-binomial/gamma-Poisson. These differences are attributable 
to the more flexible covariance structure induced by the latter, i.e., allowing for different 
covariances between the repeated number of trials. 

The values of the AIC ( = 1888.0) and the BIC ( = 1919.1) for the beta-binomial/gamma- 
Poisson model compared to the corresponding values (AIC = 1935.6 and BIC = 1974.0) 
for the multivariate beta-binomial/Poisson also suggest a better fit of the former. 

Although the results are quite similar, with the exception of the values for patients in 
the advanced stage of the disease, the beta-binomial/gamma-Poisson one is preferable to 
the multivariate beta-binomial/Poisson, both because of the modelling flexibility and the 
computational advantages mentioned before. 

As an extension for the beta-binomial/gamma-Poisson model, we could incorporate a 
parameter to relate the probabilities of success to the total attempts, as in Zhu et al. 
(2004). Another possible extension would be to consider the case where attempts could be 
done correctly, satisfactorily or incorrectly; in this case, we could generalize the model by 
considering Dirichlet-multinomial/gamma-Poisson distribution models. These extensions 
are currently under investigation. 
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Appendix 

First and second derivatives for the gamma-Poisson model 



dL((3 x ,(3 a ,(3 s 



Z' A L [L _1 n - (I p ® B" 1 )(l p ® a)] , 



dL(/3 x l3 a ,/3 s ) = KM[c _ D -i log(b)] and 
dL((3 x P a ,(3 5 ) = z , [De + D _i Mlog(b) _ B _i Lsa] 



d 2 L{(3 x ,(3 a ,(3 s ) 



Z A L[I p ® (AB" 1 )] {L[I P ® (DB" 1 )] - I A/p } Z A , 



d/3 A d/3 
d 2 L((3 x ,(3 a ,(3 s ) 



d 2 L((3 x ,(3 a ,(3 s ) 



Z A L[I P ® (DB" 2 )] [I p ® (N s - ML,)] (lp ® Z 6 



Z' Q M [C - D -1 log(B) - MF] Z c 



d 2 L((3 x ,P a ,f3 5 ) 



Z' a MB [-3 - D^B % + D" 2 log(B)] Z 5 and 



d 2 L(P x ,(3 a ,(3 s 



d(3 a d(3' s 

- = Z' 5 D {E - DQ + MD _1 B _1 L S - D~ 2 Mlog(B) - B _2 L S [N s - ML S ]} Z s 



d(3 s d(3' 5 

with L(j3 x ,(3 a ,{3 s ) presented in (24) and 



a = (ai, a 9 , ajvf )', a 9 = (5(z 5ff ) 
A = diag{a g } 



.h=l 



+ "(Zag) 



J^A(z A9h ) 



./i=i 



B = diag{b g }, b g = S(z Sg ) 

log(b) = (logih), ...,log(b g ), ...,log(b M ))' 
log(B) = diag{log(b g )} 

E L=i n sh -1 

C = (Cl, ...,Cg, ...,C M )', Cg = ^ 

C = diag{c g } 



+ 1 



u=0 



tt(Zag) +U(5(z 5g ) ! 
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u 



e = (e u ...,e g ,...,e M y, e 9 = g ^ Zag) + u6{zSg) 
E = diag{e g }, 

Z p h=1 n gh -l 

J = diag{j g }, j g = V] ^ — — ^ 



S h=l n 9& -1 



Q = diag{q g }, q g = ^ 

n = (mi, ...,n gh , ...,n Mp )', 

N s = diag j^/V' j 

L = dia5{A(z A g/ l )} 

L s = diag | g \(z Xgh ) | 

M = diag{a(z aff )} 
D = diag{8(z Sg )} 

^A = ( Z A11) • •) z A 9 /n • ) z Wp) 

— \ z ali ••■) z agi ••■) z Q Af J 

Z/ / / / \/ 

5 = { z SlT--i z Sq,--; z 8Al) 



It 



1 2 



a(z ag ) + u5(z Sg ) 
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